There are different ways to plot maps in R. The packages you use will depend on which your aim is and on the format of your data. Sometimes, geographical data can be simply in a .txt or .csv file, where you have one column for latitude and one for longitude and other variables or measurements corresponding to that location in the following columns. Yo can also have shape files, which are special files for geographical information. This files may contain points, lines or polygons. They are usually .shp files. Finally you can have raster files which are something similar to pictures or .png files because these files have pixels and (at least) a value associated with each pixel. The difference between a .png file and a shape file is that the shape-file has associated coordinates that locate this pixels on a specific surface of the world. The format is usually .geotiff or just .tiff.
In this meetup we will cover the basis of working with this formats. But this is only the starting point. There is a lot out there.
We can just use geom_polygon from the library ggplot to create our first figure. Let´s suppose we know the latitude and the longitude of all the nodes in our polygon. How can we plot it?
require(ggplot2)
funny <- data.frame(lat = c(62, 55, 48, 48, 62), long = c(20,10.5,20,10,10)) # from top right
#following the clock
ggplot(funny, aes(x = long, y = lat)) + # we specify the data
geom_polygon(fill="green") + # we plot it
geom_point(aes(x=13.25, y=50.5), color="violet", size=18) # we can also add points this was
The geometry of the polygons of the countries in the world are already loaded into R. You can access this info using the function map_data from the package ggplot2.
How do we plot only one or a subset of countries?
We need to get the data for the geometry of the countries. To do that we use the library maps and ggplot2.
countries <- c("Germany", "France")
# Get the data
some.maps <- map_data("world", region = countries) # This function retrieves the data. The data we get is a df with lat and long of each node.
head(some.maps)
## long lat group order region subregion
## 1 14.21367 53.87075 1 1 Germany Usedom
## 2 14.17217 53.87437 1 2 Germany Usedom
## 3 14.04834 53.86309 1 3 Germany Usedom
## 4 13.92578 53.87905 1 4 Germany Usedom
## 5 13.90215 53.93896 1 5 Germany Usedom
## 6 13.92168 53.99663 1 6 Germany Usedom
require(dplyr)
require(tidyr)
# Mean of lat and long to then write the labels
lab.data <- some.maps %>% group_by(region) %>% summarise(long = mean(long),
lat = mean(lat))
lab.data
## # A tibble: 2 x 3
## region long lat
## <chr> <dbl> <dbl>
## 1 France 3.23 46.2
## 2 Germany 10.4 51.2
Themes in ggplot In ggplot2 the display of non-data components is controlled by the theme system. The theme system of ggplot2 allows the manipulation of titles, labels, legends, grid lines and backgrounds. There are various build-in themes available that already have an all-around consistent style. A very common one for plotting data is theme_bw() and theme_map(). We use theme_map() which will give our plot the appearance of a map. To use theme_map() we need to load the package ggthemes().
require(ggthemes)
fra_and_ger <- ggplot(some.maps, aes(x = long, y = lat)) +
geom_polygon(aes( group = group, fill = region))+ #plot polygon, use group to group all the nodes of the same country together.
geom_text(aes(label = region), data = lab.data, size = 3, hjust = 0.5)
#using theme_map
fra_and_ger + #print labels
theme_map() +
theme(legend.position = "none")
#Using theme bw
fra_and_ger + #print labels
theme_bw() +
theme(legend.position = "none")+
xlab("Longitude") + ylab("Latitude")
ggplot2 is becoming the standard for R graphs. However, it does not handle spatial data specifically. Handling spatial objects in R relies on Spatial classes defined in the package sp or sf. ggplot2 allows the use of spatial data through the package sf. sf elements can be added as layers in a graph. The combination of ggplot2 and sf enables to create maps, using the grammar of graphics,but incorporation geographical info.
ggplot2.The function ne_countries also retrieves info about the countries. But now you do not get just nodes with coordinates but an sf object that contains much more info.
require(rnaturalearth)
require(rnaturalearthdata)
require(sf)
world <- ne_countries(scale = "medium", returnclass = "sf") # this is another function to get polygons of countries.
class(world)
## [1] "sf" "data.frame"
dim(world)
## [1] 241 64
head(world[, 1:5])
## Simple feature collection with 6 features and 5 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -70.06611 ymin: -18.01973 xmax: 74.89131 ymax: 60.40581
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## scalerank featurecla labelrank sovereignt sov_a3
## 0 3 Admin-0 country 5 Netherlands NL1
## 1 1 Admin-0 country 3 Afghanistan AFG
## 2 1 Admin-0 country 3 Angola AGO
## 3 1 Admin-0 country 6 United Kingdom GB1
## 4 1 Admin-0 country 6 Albania ALB
## 5 3 Admin-0 country 6 Finland FI1
## geometry
## 0 MULTIPOLYGON (((-69.89912 1...
## 1 MULTIPOLYGON (((74.89131 37...
## 2 MULTIPOLYGON (((14.19082 -5...
## 3 MULTIPOLYGON (((-63.00122 1...
## 4 MULTIPOLYGON (((20.06396 42...
## 5 MULTIPOLYGON (((20.61133 60...
We can plot this geographical object using ggplot2 and geometry geom_sf
plot_w1 <- ggplot(data = world) + theme_bw() + geom_sf() + xlab("Longitude") +
ylab("Latitude") + ggtitle("World map", subtitle = paste0("(", dim(world)[1],
" countries)"))
plot_w1
This really looks like a map!
plot_w1 +
geom_sf(color="blue", fill="black") +
theme(panel.background = element_rect(fill = "black"))+# Modify the theme to change the background
#Where was the funny polygon located?
geom_polygon(data=funny, aes(x = long, y = lat), color="green", fill="green") + # we plot it
# Add points using lat-lon information
geom_point(aes(x=13.25, y=50.5), color="violet", size=2)
However our nice modern-art-style polygon looks different in the map compared to the first plot. Any idea about what is happening?
You can get a feeling of how the Mercator projection distorts our worldview at link.
Google and many apps use unprojected coordinates. When the coordinates are unprojected, they are in degrees and give the position on an sphere and not in a 2D surface. But, you still need an ellipsoid and datum to make clear which reference system you are using. All this information is coded in an EPSG code. The most common coordinate reference system crs (used by Google and most apps) is EPGS: 4326. When you have lat Lon coordinates and no more info, it is likely that the EPSG is 4326, which uses the datum WGS84.
To get a taste: EPSG: 4326 Coordinate Systems Worldwide
To find the one used in your region of interest:
What is a Chorophlet Map? “A choropleth map is a thematic map in which areas are shaded or patterned in proportion to the measurement of the statistical variable being displayed on the map, such as population density or per-capita income.”(Wikipedia)
require(readr) # to use the function read_csv()
life.exp <- read_csv("data/LifeExp.csv")
The data was obtained from link, lot of cool data-sets are available for free.
life_exp <- life.exp %>%
filter(year == 2016 & sex == "Both sexes") %>% # Keep data for 2016 and for both sex
dplyr::select(country, value) %>% # Select the two columns of interest
rename(name = country, lifeExp = value) %>% # Rename columns
#We have a very common proble when using different sources, the name of the countries not allways match
# Replace "United States of America" by USA in the region column
mutate(
name = ifelse(name == "United States of America", "United States", name),
name = ifelse(name == "Russian Federation", "Russia", name)
)
# attribute names
colnames(world)
## [1] "scalerank" "featurecla" "labelrank" "sovereignt" "sov_a3"
## [6] "adm0_dif" "level" "type" "admin" "adm0_a3"
## [11] "geou_dif" "geounit" "gu_a3" "su_dif" "subunit"
## [16] "su_a3" "brk_diff" "name" "name_long" "brk_a3"
## [21] "brk_name" "brk_group" "abbrev" "postal" "formal_en"
## [26] "formal_fr" "note_adm0" "note_brk" "name_sort" "name_alt"
## [31] "mapcolor7" "mapcolor8" "mapcolor9" "mapcolor13" "pop_est"
## [36] "gdp_md_est" "pop_year" "lastcensus" "gdp_year" "economy"
## [41] "income_grp" "wikipedia" "fips_10" "iso_a2" "iso_a3"
## [46] "iso_n3" "un_a3" "wb_a2" "wb_a3" "woe_id"
## [51] "adm0_a3_is" "adm0_a3_us" "adm0_a3_un" "adm0_a3_wb" "continent"
## [56] "region_un" "subregion" "region_wb" "name_len" "long_len"
## [61] "abbrev_len" "tiny" "homepart" "geometry"
life_exp2 <- left_join(world, life_exp, by = "name")
life_exp_map <- ggplot(data = life_exp2) + geom_sf(aes(fill = lifeExp)) + scale_fill_viridis_c(option = "plasma",
trans = "sqrt") # this allows you to choose different colour scale
life_exp_map
Clearly we do not have information for all the countries regarding life expectancy and/or we have some problems with non-matching country names. When we call left_join we only keep in the table the countries that are in the left table.
We can also add points to our plot. For example we could add some capitals. Points where samples were obtained, etc. We can also control the size and color of the points using another variable. For example, we could plot sample points and the size would be number of observations. We could plot one point per city and the size number of inhabitants…and so son and so forth.
country_capitals <- read_csv("data/country-capitals.csv")
south_america_capitals <- country_capitals %>% dplyr::filter(ContinentName == "South America")
south_america_capitals$CapitalLatitude <- as.numeric(south_america_capitals$CapitalLatitude)
require(ggrepel)
life_exp_map +
geom_point(data=south_america_capitals,
aes(x=CapitalLongitude, y=CapitalLatitude),
alpha=0.5)+
geom_text_repel(data=south_america_capitals, # geom_text_repel avoidsoverlapping
aes( x=CapitalLongitude, y=CapitalLatitude, label=CapitalName),
hjust=0, vjust=0, size= 3)+
theme_bw()
require(ggspatial)
life_exp_map +
annotation_scale(location = "bl", width_hint = 0.5) + # scale
annotation_north_arrow(location = "bl", which_north = "true", #north arrow
pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering) +
coord_sf(xlim = c(-20.15, 50.12), ylim = c(35.65, -45)) #crop
Usually spatial data comes as a raster or as a polygon format. We usually have shape files with the onformation of sampling sites, neighborhoods, streets, rivers, hospitals, surveys, etc. saved as .shp files. Or we have lat-Lon data that we can transform into an sf object as we did before. We can also get raster files such as climate raster files or remote sensing data, land use, etc. Most of the times we want to plot all this info together (raster + one or more shape-files).
require(rgdal)
require(raster) #Package to work with raster files
require(sf) #Package to work with shape files
All this libraries allow us to work with spatial data using R.
rgdal is a library that that allows us to read and write geospatial data in R. The library just translates the already existing library gdal into R. The link to gdalproject is link.
raster is the library that we use to work with raster formats.
sf encodes spatial vector data. We already used it.
Let´s suppose he have different sampling points in Germany. We are interested in the temperature difference between summer and winter in this sites. Can we do a map to display this info? How would you do that?
Link to the Federal Agency for Cartography Link to open BW data Link to WoldClim link to land use data Another useful link
All the data used in this example was download from this sites or directly from R.
t_july <- raster("data/wc2.0_5m_tavg_07.tif") # loads the raster
t_december <- raster("data/wc2.0_5m_tavg_01.tif") # loads the raster
# Get temperature range
t_diff <- abs(t_july - t_december)
plot(t_diff, main = "Temperature range")
Raster calculation is relatively straight forward using R if the raster have the same resolution. If the raster come from different sources you usually have to re-shape one raster.
more info here
sites <- st_read("data/sample_sites.shp") # loads the vector data - sites
## Reading layer `sample_sites' from data source `C:\Users\ichas\Documents\RLadies\geoR\rmarkdown\data\sample_sites.shp' using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 3 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 8.081778 ymin: 47.73776 xmax: 13.30927 ymax: 53.69961
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
# get germany voundaries from world R data
germany <- world %>% dplyr::filter(name == "Germany")
plot(germany)
As you can see here, each polling has associated data. Geometry has the info to build the polygon. level, type, etc. is the info related to each polygon. If the Polygons have different info in this fields, the polygons will be displayed in different colors.
# get germany voundaries from world R data
gfs <- world %>% dplyr::filter(name %in% c("Germany", "Austria"))
plot(gfs)
We have data coming from different sources, we have to check that the coordinate reference system match (otherwise we will have lots of problems and we will be doing everything wrong)
# The function to get the coor. system is different between shapes and
# rasters
st_crs(germany) #get the projection
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
st_crs(sites) #get the projection
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
crs(t_diff) #get the projection
## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
st_crs(germany)$proj4string == st_crs(sites)$proj4string
## [1] TRUE
st_crs(germany)$proj4string == crs(t_diff, asText = T)
## [1] FALSE
The data do not have the same coord. reference system. We have to transform one.
sites_transformed <- st_transform(sites, crs = crs(t_diff, asText = T)) #transform coordinate sistem
germany_transformed <- st_transform(germany, crs = crs(t_diff, asText = T)) #transform coordinate
We get the extent (the coordinates of the extreme points) for our shape file of Germany.
germany_extent <- extent(germany_transformed)
germany_extent
## class : Extent
## xmin : 5.85752
## xmax : 15.0166
## ymin : 47.27881
## ymax : 55.05874
# crop the land use data by the extend of BW. Crop funciton from the raster
# package
crop_tdiff <- crop(t_diff, germany_extent, snap = "out")
plot(crop_tdiff)
In order to get a raster that matches exactly with the shape, we have to mask it. We could also do this step without cropping first, but that would be much more computational intensive.
mask_tdiff <- mask(crop_tdiff, germany_transformed)
plot(mask_tdiff)
tmapThe library tmap is another grate option to build thematic maps. I find it very useful to work with raster data. Useful info can be found in the tmap vignette. Another library to plot raster data is rasterVis. We will not cover it today but you can check the vignette You can plot many layers together. Each needs some reference(point, raster, text, etc.)
require(tmap)
tmap_mode("plot") # static plot or leaflet.
tm_shape(mask_tdiff) +
tm_raster (col= "layer" , style = "cont", n=10, title = "T. Difference") +
tm_shape(germany_transformed) + #add ploygon of germany
tm_borders() +
tm_shape(sites_transformed) +
tm_symbols(col = "gray", scale = .5)+ ##add points
tm_text("Site", size = 0.8, ymod=-0.6, root=1, size.lowerbound = .60, #add text
bg.color="gray", bg.alpha = .5) +
tm_layout(inner.margins = c(0.15, 0.30, 0.10, 0.1)) + # crop the extent of the map
tm_layout(legend.position = c("left","bottom")) + # add legent
tm_compass() + # add north
tm_scale_bar() #add scale
tm <- tm_shape(mask_tdiff) +
tm_raster (col= "layer" , style = "cont", n=10, title = "T. Difference") +
tm_shape(germany_transformed) + #add ploygon of germany
tm_borders() +
tm_shape(sites_transformed) +
tm_symbols(col = "gray", scale = .5)+ ##add points
tm_text("Site", size = 0.8, ymod=-0.6, root=1, size.lowerbound = .60, #add text
bg.color="gray", bg.alpha = .5) +
tm_layout(inner.margins = c(0.15, 0.30, 0.10, 0.1)) + # crop the extent of the map
tm_layout(legend.position = c("left","bottom")) + # add legent
tm_compass() + # add north
tm_scale_bar() #add scale
tmap_save(tm, filename = "t_range.png")
With tmap is really easy to do interactive maps!
require(tmap)
tmap_mode("view")
tm_shape(mask_tdiff) + tm_raster(col = "layer", style = "cont", n = 10, title = "T. Difference") +
tm_shape(germany_transformed) + tm_borders() + tm_shape(sites_transformed) +
tm_symbols(col = "gray", scale = 0.5) + tm_text("Site", size = 0.8, ymod = -0.6,
root = 1, size.lowerbound = 0.6, bg.color = "gray", bg.alpha = 0.5) + tm_layout(inner.margins = c(0.1,
0.3, 0.1, 0.1)) + tm_layout(legend.position = c("left", "bottom"))
## Reclasify values in the rasta flie. Reclasify function of raster file
## function
tmap is a grate library to make maps. Check the vignette
Get eh data for the following countries (or a subset) :Switzerland, Germany, Austria, Belgium, Netherlands, Denmark, Poland, Italy, Croatia, Slovenia, Hungary, Slovakia,Czech republic.
Plot each country with a different color.
Add the name in the center of the country.
require(ggplot2)
require(tidyr)
require(dplyr)
# Some EU Contries
some_eu_countries <- c( "Switzerland", "Germany",
"Austria", "Belgium", "Netherlands",
"Denmark", "Poland", "Italy",
"Croatia", "Slovenia", "Hungary", "Slovakia",
"Czech republic"
)
# Retrievethe map data
some_eu_maps <- map_data("world", region = some_eu_countries)
# Compute the centroid as the mean longitude and lattitude
# Used as label coordinate for country's names
region_lab_data <- some_eu_maps %>%
group_by(region) %>%
summarise(long = mean(long), lat = mean(lat))
ggplot(some_eu_maps, aes(x = long, y = lat)) +
geom_polygon(aes( group = group, fill = region)) +
geom_text(aes(label = region), data = region_lab_data, size = 3, hjust = 0.5) +
scale_fill_viridis_d() +
theme_void()+ # This is what is defining the background and the colors
theme(legend.position = "none")
require(rnaturalearth)
require(rnaturalearthdata)
require(sf)
require(ggplot2)
require(ggthemes)
require(ggspatial)
world <- ne_countries(scale = "medium", returnclass = "sf") # this is another function to get polygons of countries.
# The names of the countries are different between the sourcecs. To check
# how to spell the names:
head(unique(world$name))
## [1] "Aruba" "Afghanistan" "Angola" "Anguilla" "Albania"
## [6] "Aland"
# Czech republic is averebiated in this case, we change it:
# Some EU Contries
some_eu_countries <- c("Switzerland", "Germany", "Austria", "Belgium", "Netherlands",
"Denmark", "Poland", "Italy", "Croatia", "Slovenia", "Hungary", "Slovakia",
"Czech Rep.")
some_eu_sf <- world %>% dplyr::filter(name %in% some_eu_countries)
# This is a possible solution
ggplot(data = some_eu_sf) + geom_sf(aes(fill = name)) + coord_sf(xlim = c(0,
25.12), ylim = c(35.65, 60)) + theme_bw()
# To get the names cetered
eu_points <- st_centroid(some_eu_sf)
eu_points <- cbind(eu_points, st_coordinates(st_centroid(eu_points$geometry)))
ggplot(data = some_eu_sf) + geom_sf(aes(fill = name)) + coord_sf(xlim = c(0,
25.12), ylim = c(35.65, 60)) + theme_map() + theme_bw() + geom_text(data = eu_points,
aes(x = X, y = Y, label = name), fontface = "bold", check_overlap = FALSE,
size = 3, color = "darkred") + theme(legend.position = "none") + xlab("Longitude") +
ylab("Latitude") + scale_fill_viridis_d() + annotation_scale(location = "bl",
width_hint = 0.3) + annotation_north_arrow(location = "bl", which_north = "true",
pad_x = unit(0.15, "in"), pad_y = unit(0.3, "in"), style = north_arrow_fancy_orienteering)
require(tmap)
new_bb = c(-10, 35, 35, 60)
names(new_bb) = c("xmin", "ymin", "xmax", "ymax")
attr(new_bb, "class") = "bbox"
attr(st_geometry(some_eu_sf), "bbox") = new_bb
tmap_mode("plot")
tm_shape(some_eu_sf) + tm_borders() + tm_polygons("name") + tm_shape(eu_points) +
tm_symbols(col = "gray", scale = 0.5) + tm_text("name", size = 0.8, ymod = -0.6,
root = 1, size.lowerbound = 0.6, bg.alpha = 0.5) + tm_compass() + tm_scale_bar() +
tm_legend(show = FALSE)
This data frame has information about all the R Ladies chapters around the world! It is quite interesting. There are many questions we can answer. But first we have to take a look at the data.
require(readr)
rladies <- read_csv("data/RLadiesChapters.csv")
head(rladies)
## # A tibble: 6 x 8
## X1 chapter date members lat lon city country
## <int> <chr> <chr> <int> <dbl> <dbl> <chr> <chr>
## 1 1 R-Ladies Ba~ 22/10/201~ 389 41.4 2.17 Barcelona ES
## 2 2 R-Ladies Us~ 09/05/201~ 23 -54.8 -68.3 Ushuaia AR
## 3 3 R-Ladies Bi~ 27/02/201~ 38 43.2 -2.93 Bilbao ES
## 4 4 R-Ladies Ba~ 10/05/201~ 114 -41.1 -71.3 San Carlos ~ AR
## 5 5 R-Ladies Me~ 02/09/201~ 1213 -37.8 145. Melbourne AU
## 6 6 R-Ladies No~ 24/05/201~ 46 45.2 19.8 Novi Sad RS
# date stands for the time when the chapter was created. We will change the
# name of this column
rladies$created_at <- as.character(rladies$date)
Which would be a nice plot to show this data? Would you use a clorophlet map? Or Would you add points? There is no right or wrong! Let´s try:
Where are the chapters of R Ladies? Where are the bigger/smaller chapters?
require(ggplot2)
require(ggthemes)
# Here we do the same map we did at the begining of the Meetup
world <- ggplot() + borders("world", colour = "gray85", fill = "gray80") + theme_map()
map <- world + # We create points
geom_point(aes(x = lon, y = lat, size = members), data = rladies, colour = "purple",
alpha = 0.5) + scale_size_continuous(range = c(1, 8), breaks = c(250, 500,
750, 1000)) + labs(size = "Followers") + ggtitle("R-Ladies chapters around the world")
map
In which countries are there more R Ladies chapters? A clorophlet map.
require(dplyr)
require(tidyr)
# We summarise the number of chapters per country
country.chapters <- rladies %>% group_by(country) %>% summarise(n.chapters = n())
# We have to add this data to the world_map data that has the poligons
world_map <- map_data("world")
However we have a problem here…we have to join the data-sets but the country names are codded in a different way. This kind of problems are really common when we work with countries. I just found a txt with country name and code online here
countries <- read.csv("data/countries.txt")
colnames(countries)[2] <- "country" #use the same name as in rladies df
rladies.country <- left_join(country.chapters, countries, by = "country") # join rladies df with the names of the countries
colnames(rladies.country)[3] <- "region"
map.rladies <- full_join(rladies.country, world_map, by = "region") # join the new rlades df with the map df
ggplot(map.rladies, aes(long, lat, group = group)) + geom_polygon(aes(fill = n.chapters),
color = "white") + scale_fill_viridis_c(option = "C", tran = "log", breaks = c(1,
3, 10, 40)) + theme_bw() + ggtitle("R-Lchapters per country")
The tif file g100_06.tif contains land unse information for Europ. The file AX_Gebiet_Kreis.shp is a shape file of Baden Württemberg. Visualize the land use for BW. Do not forget to check crs. I have been trying, but I could not make the legend appear. Info for the legend in is the file clc_legend.csv. May be you find out!
# Load data
require(raster)
require(sf)
corine <- raster("data/g100_06.tif") # loads the raster
shape <- st_read("data/AX_Gebiet_Kreis.shp") # loads the vector data
## Reading layer `AX_Gebiet_Kreis' from data source `C:\Users\ichas\Documents\RLadies\geoR\rmarkdown\data\AX_Gebiet_Kreis.shp' using driver `ESRI Shapefile'
## Simple feature collection with 44 features and 2 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 388337.8 ymin: 5265159 xmax: 610069.4 ymax: 5515632
## epsg (SRID): NA
## proj4string: +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs
# Inspect and plot the data Different functions to see the info that is
# encoded in our shape
dplyr::glimpse(shape)
## Observations: 44
## Variables: 3
## $ Name <fct> Lörrach, Waldshut, Bodenseekreis, Ravensburg, Breisg...
## $ Schlüssel <fct> 08336, 08337, 08435, 08436, 08315, 08326, 08327, 083...
## $ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((418041.5 53..., MULTIPO...
colnames(shape)
## [1] "Name" "Schlüssel" "geometry"
summary(shape)
## Name Schlüssel geometry
## Heilbronn : 2 08111 : 1 MULTIPOLYGON :44
## Karlsruhe : 2 08115 : 1 epsg:NA : 0
## Alb-Donau-Kreis: 1 08116 : 1 +proj=utm ...: 0
## Baden-Baden : 1 08117 : 1
## Biberach : 1 08118 : 1
## Böblingen : 1 08119 : 1
## (Other) :36 (Other):38
plot(shape)
# Plot the raster
plot(corine)
# Check coordinate reference system The function to get the coor. system is
# different between shapes and rasters
st_crs(shape) #get the projection
## Coordinate Reference System:
## No EPSG code
## proj4string: "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"
crs(corine) #get the projection
## CRS arguments:
## +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
st_crs(shape)$proj4string == crs(corine, asText = T)
## [1] FALSE
shp.Transformed <- st_transform(shape, crs = crs(corine, asText = T)) #transform coordinate sistem
plot(st_geometry(shp.Transformed)) # plots only the geometry and not the associated values
st_crs(shp.Transformed)$proj4string == crs(corine, asText = T) #check again if they are the same coordinate sistem
## [1] TRUE
# For extraction the shape file has to be of class 'Spatial Polygons Data
# Frame'
shp.Transformed.sp <- as(shp.Transformed, "Spatial")
class(shp.Transformed.sp)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
We get the extent (the coordinates of the extreme points) for our shape file of Baden Württemberg.
# get extent of Baden-Württemberg
corine.shp.extent <- extent(shp.Transformed.sp)
corine.shp.extent
## class : Extent
## xmin : 4134156
## xmax : 4357500
## ymin : 2715907
## ymax : 2964364
# We crop the land use data by the extend of BW. Crop funciton from the
# raster package
corine.shp <- crop(corine, corine.shp.extent, snap = "out")
plot(corine.shp)
# mask BW - takes about 10 minutes
corine.shp <- mask(corine.shp, shp.Transformed.sp)
# 7) save corine.bav as a new file
writeRaster(corine.shp, filename = "data/corineBW.tif", format = "GTiff", overwrite = TRUE)
# loads it again
corine.shp <- raster("data/corineBW.tif")
plot(corine.shp, main = "Land Use")
scalebar(10, xy = click(), type = "line")
require(tmap)
# Visualize 1) import legend from a csv with readr package
legend <- readr::read_csv2("data/clc_legend.csv")
# note a special feature: attribute progress
`?`(readr::read_csv2)
# how does the legend look like?
knitr::kable(head(legend))
| GRID_CODE | CLC_CODE | LABEL1 | LABEL2 | LABEL3 | RGB |
|---|---|---|---|---|---|
| 1 | 111 | Artificial surfaces | Urban fabric | Continuous urban fabric | 230-000-077 |
| 2 | 112 | Artificial surfaces | Urban fabric | Discontinuous urban fabric | 255-000-000 |
| 3 | 121 | Artificial surfaces | Industrial, commercial and transport units | Industrial or commercial units | 204-077-242 |
| 4 | 122 | Artificial surfaces | Industrial, commercial and transport units | Road and rail networks and associated land | 204-000-000 |
| 5 | 123 | Artificial surfaces | Industrial, commercial and transport units | Port areas | 230-204-204 |
| 6 | 124 | Artificial surfaces | Industrial, commercial and transport units | Airports | 230-204-230 |
## 4) plot map using tmap package
tmap_mode("plot")
tm_shape(corine.shp) + tm_raster(style = "cat", title = "Categories") + tm_layout(legend.position = c("left",
"bottom")) + tm_shape(shape) + tm_borders()
Here we will start justo from a data frame. The df has info of bird observations of one species (Pyrocephalus rubinus) in South and Central America. This bird is migratory, we want to visualize the location of the observations in different moments of the year. How would you do it?
The df is quite big, so to speed up I woud suggest to use recent observatins, e.g. after 2016.
require(lubridate)
require(ggplot2)
require(sf)
require(ggspatial)
require(dplyr)
require(tidyr)
require(rnaturalearth)
require(rnaturalearthdata)
load("data/Pyrocephalus-rubinus.RData")
Py_rubunis <- birds_df %>% mutate(m = month(observation_date), y = year(observation_date)) %>%
dplyr::select(m, y, latitude, longitude, species_observed) %>% dplyr::filter(!is.na(latitude)) %>%
dplyr::filter(!is.na(longitude)) %>% dplyr::filter(m %in% c(1, 7)) %>% dplyr::filter(species_observed ==
TRUE) %>% dplyr::filter(y > 2016)
# We get
crs <- "+proj=longlat +datum=WGS84 +no_defs" #EPSG 4326
Py_rubunis_geo <- st_as_sf(Py_rubunis, coords = c("longitude", "latitude"),
crs = crs)
world <- ne_countries(scale = "medium", returnclass = "sf")
ggplot(data = world) + geom_sf() + theme_bw() + geom_sf(data = Py_rubunis_geo,
aes(color = as.factor(m))) + xlab("Longitude") + ylab("Latitude") + coord_sf(xlim = c(-120.15,
-25), ylim = c(40, -75)) + ggtitle("Observations of P. rubinus in January and July")
load("data/Pyrocephalus-rubinus.RData")
require(tmap)
require(sf)
require(rnaturalearth)
require(rnaturalearthdata)
require(dplyr)
require(tidyr)
require(lubridate)
world <- ne_countries(scale = "medium", returnclass = "sf") # this is another function to get polygons of countries.
countries <- c("Argentina", "Bolivia", "Brazil", "Chile", "Colombia", "Ecuador",
"Guyana", "Paraguay", "Peru", "Suriname", "Uruguay", "Venezuela", "Mexico",
"Belize", "Guatemala", "El Salvador", "Honduras", "Nicaragua", "Costa Rica",
"Panama")
sc_america <- world %>% dplyr::filter(name %in% countries)
Py_rubunis <- birds_df %>% mutate(m = month(observation_date), y = year(observation_date)) %>%
dplyr::select(m, y, latitude, longitude, species_observed) %>% dplyr::filter(!is.na(latitude)) %>%
dplyr::filter(!is.na(longitude)) %>% dplyr::filter(species_observed == TRUE) %>%
dplyr::filter(y > 2016)
crs <- "+proj=longlat +datum=WGS84 +no_defs" #EPSG 4326
Py_rubunis_geo <- st_as_sf(Py_rubunis, coords = c("longitude", "latitude"),
crs = crs)
Py_rubunis_1 <- Py_rubunis %>% dplyr::filter(m %in% c(1))
Py_rubunis_4 <- Py_rubunis %>% dplyr::filter(m %in% c(4))
Py_rubunis_7 <- Py_rubunis %>% dplyr::filter(m %in% c(7))
Py_rubunis_10 <- Py_rubunis %>% dplyr::filter(m %in% c(10))
Py_rubunis_geo_1 <- st_as_sf(Py_rubunis_1, coords = c("longitude", "latitude"),
crs = crs)
Py_rubunis_geo_4 <- st_as_sf(Py_rubunis_4, coords = c("longitude", "latitude"),
crs = crs)
Py_rubunis_geo_7 <- st_as_sf(Py_rubunis_7, coords = c("longitude", "latitude"),
crs = crs)
Py_rubunis_geo_10 <- st_as_sf(Py_rubunis_10, coords = c("longitude", "latitude"),
crs = crs)
tmap_mode("view")
tm_shape(sc_america) + tm_borders() + tm_polygons() + tm_shape(Py_rubunis_geo_1) +
tm_symbols(col = "blue", scale = 0.5) + tm_shape(Py_rubunis_geo_4) + tm_symbols(col = "green",
scale = 0.5) + tm_shape(Py_rubunis_geo_7) + tm_symbols(col = "orange", scale = 0.5) +
tm_shape(Py_rubunis_geo_10) + tm_symbols(col = "red", scale = 0.5)
# tm_compass() + tm_scale_bar() + tm_legend(show = FALSE)